MODULE md1_basis_set

   USE md1_global_types
   USE md1_global_consts
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: md1_init_basis_set,         &
             md1_free_basis_set,         &
             md1_get_basis_set,          &
             md1_get_basis_set_extent

   TYPE(md1_basis), TARGET, SAVE :: basis
   LOGICAL, SAVE :: md1_basis_set_init = .FALSE.

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE md1_init_basis_set(filename)

      IMPLICIT NONE
      CHARACTER(*), INTENT(IN) :: filename

      CALL init_Ltuv_indices
      CALL init_raw_data(filename)
      CALL init_AO_indices
      md1_basis_set_init = .TRUE.

   END SUBROUTINE md1_init_basis_set

!-------------------------------------------------------------------------------

   SUBROUTINE md1_get_basis_set(basis_ptr)

      IMPLICIT NONE
      TYPE(md1_basis), POINTER :: basis_ptr

      IF (md1_basis_set_init) THEN
         basis_ptr => basis
      ELSE
         STOP 'md1_basis_set not initialized!'
      END IF

   END SUBROUTINE md1_get_basis_set

!-------------------------------------------------------------------------------

   SUBROUTINE md1_get_basis_set_extent(thresh,AO_extent)

      USE md1_shell_pairs, ONLY: md1_get_shell_pairs
      USE md1_utils,       ONLY: md1_get_prim_batch

      IMPLICIT NONE
      REAL(REALK),         INTENT(IN)  :: thresh
      TYPE(md1_real_hilo), INTENT(OUT) :: AO_extent(3)

      TYPE(md1_sh_pairs), POINTER :: sh_pairs(:)
      TYPE(md1_prim_batch) :: batch(basis%Cmax2)
      INTEGER(INTK) :: i,k, xyz, np
      REAL(REALK) :: lo, hi, r
      
      ! Get list of non-vanishing shell pairs
      CALL md1_get_shell_pairs(basis,sh_pairs)

      ! Find min coordinates at which all basis functions decay to threshold
      lo =  1d10
      hi = -1d10
      DO k = 1, SIZE(sh_pairs)
         CALL md1_get_prim_batch(basis,sh_pairs(k),batch,np)
         DO xyz = 1, 3
            DO i = 1, np
               r  = SQRT(-LOG(thresh) / batch(i)%exp)
               lo = MIN(lo, batch(i)%P(xyz) - r)
               hi = MAX(hi, batch(i)%P(xyz) + r)
            END DO
            AO_extent(xyz)%lo = lo
            AO_extent(xyz)%hi = hi
         END DO
      END DO

   END SUBROUTINE md1_get_basis_set_extent

!-------------------------------------------------------------------------------

   SUBROUTINE init_raw_data(filename)

      IMPLICIT NONE
      CHARACTER(*), INTENT(IN) :: filename
      INTEGER(INTK) :: i,j,k, L, nprim, nbas

      OPEN (UNIT=IO2, FILE=filename, STATUS='OLD',      &
            ACTION='READ', FORM='FORMATTED')
      REWIND (IO2)

      READ (IO2,'(I8)') basis%natom
      READ (IO2,'(I8)') basis%nshell
      READ (IO2,'(I8)') nbas
      READ (IO2,'(I8)') nprim

      ALLOCATE(basis%xyz(3,basis%natom))
      ALLOCATE(basis%atom(basis%nshell))
      ALLOCATE(basis%L(basis%nshell))
      ALLOCATE(basis%exp_id(basis%nshell))
      ALLOCATE(basis%nprim(basis%nshell))
      ALLOCATE(basis%exp(nprim))
      ALLOCATE(basis%C(nprim))

      DO i = 1, basis%natom
         READ (IO2,'(F25.15)') basis%xyz(1,i) 
         READ (IO2,'(F25.15)') basis%xyz(2,i) 
         READ (IO2,'(F25.15)') basis%xyz(3,i) 
      END DO

      k = 0
      DO i = 1, basis%nshell
         READ (IO2,'(I8)') basis%atom(i)   
         READ (IO2,'(I8)') basis%L(i)   
         READ (IO2,'(I8)') basis%exp_id(i)  
         READ (IO2,'(I8)') basis%nprim(i)   
         DO j = 1, basis%nprim(i)
            k = k+1
            READ (IO2,'(F25.15)') basis%exp(k)   
            READ (IO2,'(F25.15)') basis%C(k)   
!DIMITRIJ
            WRITE(IO1,*) k, basis%exp(k)   , basis%C(k)   

         END DO
      END DO

      basis%nbas = 0
      basis%LMAX = 0
      basis%Cmax2 = 0
!FIXME change default primary basis to be spherical
! (only integrals done internally in Cartesian basis)
      DO i = 1, basis%nshell
         L = basis%L(i)
         basis%nbas = basis%nbas + basis%LM(L)%car_nlm
         basis%LMAX = MAX(basis%LMAX, L)
         basis%Cmax2 = MAX(basis%Cmax2, basis%nprim(i))
      END DO
      basis%Cmax2 = basis%Cmax2 **2
      IF (basis%nbas /= nbas) THEN
         WRITE(IO1,*) "Inconsistent number of basis funcs in basis.qcim"
         WRITE(IO1,*) "file and qcim function argument... Probably means"
         WRITE(IO1,*) "you are not using Cartesian GTOs."
         WRITE(IO1,*) "basis%nbas, nbas", basis%nbas, nbas
         STOP 'md1: init_raw_data'
      END IF

      WRITE(IO1,*), 'basis%natom  =', basis%natom
      WRITE(IO1,*), 'basis%nshell =', basis%nshell
      WRITE(IO1,*), 'basis%nbas   =', basis%nbas
      WRITE(IO1,*), 'basis%LMAX   =', basis%LMAX
      WRITE(IO1,*), 'basis%Cmax2  =', basis%Cmax2

      CLOSE (UNIT=IO2, STATUS='KEEP')

   END SUBROUTINE init_raw_data

!-------------------------------------------------------------------------------

   SUBROUTINE init_Ltuv_indices

      IMPLICIT NONE
      INTEGER(INTK) :: L,M, Lt, Lu, Lv, Maxtuv

!      Maxtuv = basis%LMAX * 2
      Maxtuv = 20  !fixme  ! depends on multipole LMAX -- actually NOT just 
                           ! a basis set issue!
 
      ALLOCATE(basis%LM(0:Maxtuv))
      DO L = 0, Maxtuv
         M = (L+1) * (L+2) / 2
         ALLOCATE(basis%LM(L)%tuv(3,M))
         basis%LM(L)%car_nlm = M
         basis%LM(L)%sph_nlm = 2*L + 1
      END DO

      ! Define order in which Cartesian components are indexed
      !              Lt   Lu   Lv
      !      L     x^   y^   z^    Ituv
      !-----------------------------------
      !      0     0    0    0      1
      !-----------------------------------
      !      1     1    0    0      2
      !            0    1    0      3
      !            0    0    1      4
      !-----------------------------------
      !      2     2    0    0      5
      !            1    1    0      6
      !            1    0    1      7
      !            0    2    0      8
      !            0    1    1      9
      !            0    0    2     10
      !-----------------------------------
      !      3     3    0    0     11
             
      DO L = 0, Maxtuv
         M = 0
         DO Lt = Maxtuv, 0, -1
            DO Lu = Maxtuv, 0, -1
               DO Lv = Maxtuv, 0, -1
                  IF (Lt+Lu+Lv == L) THEN
                     M = M + 1
                     basis%LM(L)%tuv(1,M) = Lt
                     basis%LM(L)%tuv(2,M) = Lu
                     basis%LM(L)%tuv(3,M) = Lv
                  END IF
               END DO
            END DO
         END DO
      END DO

   END SUBROUTINE init_Ltuv_indices

!-------------------------------------------------------------------------------

   SUBROUTINE init_AO_indices
   ! Define starting indices of shells of contracted Cartesian AOs
   ! e.g. AO Fock matrix p or q indices

      IMPLICIT NONE
      INTEGER(INTK) :: i

      ALLOCATE(basis%car_AO(basis%nshell))
      ALLOCATE(basis%sph_AO(basis%nshell))

      basis%car_AO(1) = 1
      basis%sph_AO(1) = 1
      DO i = 2, basis%nshell
         basis%car_AO(i) = basis%car_AO(i-1) + basis%LM(basis%L(i-1))%car_nlm
         basis%sph_AO(i) = basis%sph_AO(i-1) + basis%LM(basis%L(i-1))%sph_nlm
      END DO


   END SUBROUTINE init_AO_indices

!-------------------------------------------------------------------------------

   SUBROUTINE md1_free_basis_set

      USE md1_shell_pairs, ONLY: md1_free_shell_pairs
      IMPLICIT NONE
      INTEGER(INTK) :: i

      DEALLOCATE(basis%xyz)
      DEALLOCATE(basis%atom)
      DEALLOCATE(basis%L)
      DEALLOCATE(basis%exp_id)
      DEALLOCATE(basis%nprim)
      DEALLOCATE(basis%car_AO)
      DEALLOCATE(basis%sph_AO)
      DEALLOCATE(basis%exp)
      DEALLOCATE(basis%C)
      DO i = 1, SIZE(basis%LM)
         DEALLOCATE(basis%LM(i-1)%tuv)
      END DO
      DEALLOCATE(basis%LM)

      CALL md1_free_shell_pairs
      md1_basis_set_init = .FALSE.

   END SUBROUTINE md1_free_basis_set

!-------------------------------------------------------------------------------

END MODULE md1_basis_set
